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NONLINEAR DYNAMICS OF ATTRACTIVE MAGNETIC BEARINGS 1 
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The nonlinear dynamics of a ferromagnetic shaft suspended by the force of attraction of 1, 2, or 4 independent 
electromagnets is presented. Each model includes a state variable feedback controller which has been designed using 
the pole placement method. The constitutive relationships for the magnets are derived analytically from magnetic 
circuit theory, and the effects of induced eddy currents due to the rotation of the journal are included using Maxwell’s 
field relations. 

A rotor suspended by four electro-magnets with closed loop feedback is shown to have nine equilibrium points 
within the bearing clearance space. As the rotor spin speed increases, the system is shown to pass through a Hopf 
bifurcation (a flutter instability). Using center manifold theory, this bifurcation can be shown to be of the subcritical 
type, indicating an unstable limit cycle below the critical speed. The bearing is very sensitive to initial conditions, and 
the equilibrium position is easily upset by transient excitation. The results are confirmed by numerical simulation. 

INTRODUCTION 

The suspension of a rotating shaft in a magnetic field without mechanical contact or lubrication is an old idea. 
Many studies concerning the feasibility of electromagnetic levitation in various applications like magnetic bearings, 
high speed ground transportation, electromagnetic dampers for vibration control, etc. can be found in the literature 
[refs. 1,2,3]. Now, with advanced electronic control, high speed magnetic bearings are a commercial reality, being 
used in grinding and polishing machinery, vacuum pumps, compressors, turbines, generators, and centrifuges [refs. 
4,5]. Electromagnetic dampers have been shown to be capable of effectively eliminating vibration at the critical speed 
associated with the first bending mode of shafts [refs. 3,6]. In addition, the dampers have been seen to suppress the 
system instability associated with the fluid film bearings. 

This paper only considers attraction systems under active control. Passive systems using permanent magnets 
in repulsion are usually incapable of generating sufficient load carrying capacities. However, the constitutive model 
presented here is applicable to electromagnets operating at sufficient speed that the dead weight is carried by eddy 
currents generating lift from the lower magnet. 

The control system studied varies the voltage to each magnet in response to the motion of the mass. Such a 
suspension system with multiple magnets is a multiple input - multiple output system. A state variable feedback 
controller, designed on the basis of pole placement technique, is used to stabilize the system and meet the design 
specifications. Unlike in a single input - single output system, the multivariable problem has a much richer structure 
and has many gain matrices yielding the same pole placement [ref. 7]. The direct method of choosing an arbitrary 
vector and reducing the multi-input - multi - output system into a single input - single output system is used to 
obtain the gain matrix. 

Research and development activity on passive, active, and combination magnetic bearing systems spans over 
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150 years. Earnshaw in 1839 demonstrated [ref. 8] mathematically that it is impossible for a pole placed in a static 
field of force to have a stable equilibrium position when an inverse square law relates force and distance. Braunbek in 
1939 carried out [ref. 8] a similar mathematical analysis specifically for suspensions in unvarying magnetic fields and 
deduced that levitation is impossible in such fields when all the materials have relative permeabilities (jn r ) greater 
than or equal to 1, but possible if materials with Hr < 1 can be introduced. It follows that it is impossible to achieve 
levitation in static magnetic fields, that is, using permanent magnets or fixed current dc electromagnets, unless part 
of the system consists of either diamagnetic material (jj. r < 1) or a superconductor (ji = 0). 

Examples of hybrid passive and active systems appeared as early as 1950 in which photoelectric positional 
feedback was used. Stiffness and damping were minimal because of the lack of sophisticated control components. 
Early experiments with fully active systems (1957) were largely thwarted by the then high cost and size of control 
system components. In time, it became clear that hybrid systems with simple electronics capable of carrying industrial 
loads depended heavily on complex mechanical dampers. Meanwhile the cost and size of sophisticated electronics 
for purely active systems continued to come down as performance rose. Since active magnetic bearings provide both 
damping and support, the choice became clear. Subsequent efforts have concentrated on active magnetic bearings. 

Several reviews of electromagnetic levitation are available in the literature [ref. 8,9], Two electromagnetic 
levitation methods have met with success: direct, position feedback control techniques; and ac modulated or indirect 
feedback methods. In the latter, the magnet inductance is part of a tuned circuit whose natural frequency depends 
on the gap between the suspended mass and magnet. This method has been used to suspend gyroscopic devices for 
inertial sensors. It suffers from high eddy current losses and a small range of stable air gaps. Analysis of ac tuned 
circuit methods may be found in Kaplan [ref. 10]. 

For heavily loaded bearings, direct feedback methods have to be used. However, dynamics has not seen as much 
attention as it has in journal bearings. Most of the available literature deals with empirical ideas and concentrates 
on reliability of the bearing, reducing the size, weight and complexity of the devices [ref. 11]. 

Insight can be gained by studying a ferromagnetic suspension (suspended mass is not moving). Moon [ref. 12] 
and Woodson and Melcher [ref. 13] are good starting points. 

When the suspended rotor is spinning under the magnets and if the rotor is not laminated, eddy currents will be 
induced in the material. The induced eddy currents create two kinds of force on the rotor: drag force which leads to 
additional power dissipation and coupling of motions of the rotor in two perpendicular directions, and repulsive force 
which tends to counter balance the attractive force. Expressions for the drag and repulsive forces can be derived 
by studying the effect of material motion on the diffusion of magnetic fields [ref. 13], starting from Maxwell’s field 
relations. Moon [ref. 12] has shown that in certain magnetic levitation configurations eddy currents can produce a 
positive or negative damping force depending on the speed. 

Although a good analytical model is not available for eddy currents due to shaft rotation, a number of authors 
calculate the eddy current effects in other geometrical configurations using a hypothetical simplified model and finite 
element methods [refs. 14,15], Several studies on linear induction motors are available [refs. 1,16] which can be 
extended to magnetic bearings by making several assumptions and manipulations. In the above studies, the authors 
calculate the induced forces by cross multiplying the current density vector, which does not take care of the attractive 
force when the moving material is ferromagnetic. In this paper , the forces are calculated using the Maxwell’s stress 
tensor approach which in one calculation gives all the forces involved. 

Matsumara [ref. 17] has derived the fundamental equations for a horizontal shaft magnetic bearing taking into 
account the rolling, pitching, and yawing of the rotor. In deriving the equations of motion, he assumes that the rotor 
consists of a laminated core and consequently no eddy currents are generated in the material. He has proposed an 
integral type control system which stabilizes the system without steady state shaft position error. 

Hebbale [ref. 18] has studied the nonlinear dynamics in terms of equilibrium points, transient response, on-set 
of instability, limit cycle size, and forced response. The material which follows was taken from [ref. 18], 


NOMENCLATURE 

a Distance of magnet pole corner (near) from center line 
A Area of magnet pole face, Linearized state matrix 
b Distance of magnet pole corner (far) from center line 
B Control input matrix, a Magnetic Flux density 

Bo Magnetic flux density under the magnet pole 

Ci Feedback gains of controller 

D Diameter of rotor 
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e Unit vectors 

E{t) Control voltage 

Ei Total voltage to magnet 

E 0i Initial voltage to magnet 

f,F Nonlinear functions 

F m ag Magnetic force 

g Acceleration due to gravity 

H Magnetic field intensity 

I i Current in coil 

j V=T 

k Wave number 

K Feedback gain matrix 

L Inductance of circuit 

Li Flux path length 

M Mass of rotor 

n\,nz Real and imaginary parts of 7 

N Number of turns in coil 

q Complex wave number 

R Resistance of coil 

si Nondimensional velocity 

Time 

V Velocity vector 

V Velocity (magnitude) 

Wd Field energy density 

Wf Field energy stored 

x shaft displacement 

X State vector 

y shaft displacement 

Z Linearized state vector 

a System constant = 1 ioN 2 A 

P System constant = L]//i rl +L 2 /fJ-r2 

y Complex number = coshjfcA +fi Q q/^ik sinhitA 

T Magnetomotive force 

S variation from equilibrium 

A Equilibrium gap length 

9 Magnet pole angle with respect to vertical 

X Flux linkage 

lx,H r Permeability 

a Electrical conductivity of rotor material 

£ Closed surface, Summation sign 

C> Magnetic flux 

£2 Spin speed of rotor 


FOUR MAGNET BEARING 

The shaft is suspended within four magnets, as shown in Fig. 1. The shaft has high permeability and high 
conductivity and is not laminated. (Lamination would serve to inhibit eddy currents.) All of the four electromagnets 
are identical, modelled as a coil of N turns on a laminated core of high permeability. The coil has resistance R and 
an initial voltage £ 0 applied. The dead weight is suspended by the difference in forces exerted by top and bottom 
magnets. All magnets carry a steady state voltage, adjusted so that at zero speed, the gap lengths under the magnets 
are equal. 

This fully energized configuration was chosen so that the magnets could generate the effect of repulsion by 
decreasing the attraction. The horizontal magnets generate equal force when the shaft is centered and the control 
system is commanding equal voltages. The shaft displacement is measured by two coordinates {x,y) as shown in 


399 


Fig. 1, measured from the center of the clearance space in the horizontal and vertical directions respectively. When 
the rotor is not spinning, the center of the rotor is at the origin and the gap lengths are all equal to A . 

If the shaft is not laminated, motion of the conducting shaft through the supporting magnetic fields will generate 
induced eddy currents. Eddy current effects cause both loss of lift and a drag force in the perpendicular direction. 
It is assumed that the rotor is spinning in the clockwise direction so that the drag force due to the top magnet acts 
in the negative x direction, that due to the right magnet acts in the y direction, and so on. 

Several geometrical assumptions are made. The rotor always remains perfectly aligned within the bearing (no 
tilting). Under small displacements the surfaces of the rotor and the magnet pole faces are assumed to remain 
parallel. Since the individual poles are located at angles ±9 relative to the vertical axis in the case of vertical 
magnets or relative to the horizontal axis in the case of horizontal magnets, it is assumed that when the rotor moves 
vertically a distance y, the change in gap length for the vertical magnets is ycos9. Similarly, when the rotor moves 
horizontally by a distance x, the change in gap length for the horizontal magnets is xcosd. Any other translational 
motion of the rotor can be written as a superposition of the motions in x and y directions. The effect of unequal 
gap lengths under a vertical magnet caused by a rotor motion in the horizontal direction or vice versa is neglected 
because the total gap length under that magnet remains constant. 

The control system studied varies the voltage to each magnet in response to the motion of the mass. Such a 
suspension system with multiple magnets is a multiple input - multiple output system. A state variable feedback 
controller, designed on the basis of pole placement technique, is used to stabilize the system and meet the design 
specifications. Unlike in a single input - single output system, the multivariable problem has a much richer structure 
and has many gain matrices yielding the same pole placement [10]. The direct method of choosing an arbitrary 
vector and reducing the multi-input - multi - output system into a single input - single output system is used to 
obtain the gain matrix. 

EQUATIONS OF MOTION 

The differential equations of motion describing the response of the rotor are given by 

My— ^^F mag ■ n y — Mg 

^ „ ( 1 ) 
Mx "tag ' n x 

where ^F mag are the vectorially combined magnetic forces from all the magnets. Expressions for resolutes of F mag 
have been derived in Appendix B for a single magnet. The reader is cautioned that thex.y coordinates in Appendix B 
are local tangential and normal directions for each individual magnet and do not correspond the the x,y coordinates 
for motion of the shaft. F mag includes the steady state attraction forces and eddy current forces (both repulsion and 
drag). Only spin velocity is assumed to generate eddy currents. Motion of the shaft ( x,y ) doesn’t generate eddy 
currents. 

The remaining differential equations are obtained from a voltage balance in each of the four electromagnetic 
circuits: 


E i = I i R i +N i ti i= 1,2, 3, 4 (2) 

where subscripts refer to the top, right, bottom, and left magnets respectively. An expression for the magnetic flux 
<p has been derived in Appendix B. 

Expressions for the magnetic forces are given in Appendices A and B. If eddy currents are neglected, then there 
are no terms in the differential equations which depend on speed. In general 

Fmag =/(x,y,/ 1 ,/ 2 ,/ 3 ,/ 4 ,0) (3) 

and 

<Pi =f(x,i,y,y,i l ,f 2 ,I 3 ,l\,n) (4) 

The general procedure is to 

(a) Determine £o« such that the equilibrium location at zero speed is centered. 

(b) Design a state variable controller by pole placement. 
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(c) Determine all steady state equilibrium locations as a function of speed. 

(d) Linearize about equilibrium locations to determine stability as a function of speed. 

Choosing eight state variables y,y,.*,.x,/i,/ 2 ,/ 3 ,and/ 4 , Eqs. (1) and (2) can be put into a standard first order 

form 


X = / (X, E, Q) 


(5) 

E = Eo + SE 


(6) 

SE = KX 


(7) 

(SE i 8E 2 8E 3 

8E 4 ) t 

(8) 

y x x I\ 1 2 

h h) T 

(9) 


where X is a 8-vector containing the state variables. 


SINGLE MAGNET CASE 

The single magnet suspension serves as a paradigm for the more complicated models which follow. Consider 
only the top magnet, with x constrained to be zero, and with no rotation ( £1= 0.0). The approximation for F mag 
from Appendix A can be used, giving F mag oc (/ 2 , ^*). Furthermore, F mag is always positive. 

This system is well known to be unstable without feedback. A complete state variable feedback would be 


{ * -x 0 "I 

* j (10) 

Fig. 2 shows how the static magnetic force is affected by the choice of C\, Obviously, there is some minimum 
value of Ci to produce a positive slope at xq = 0.001. High values of Ci produce a stiffer system but with a penalty. 
A second ( unstable ) equilibrium point exists. Increasing Ci moves this point closer to the stable operating location. 
The safe operating domain in the state space is therefore decreased. 

In addition, Fig. 2.5 shows the extreme nonlinearity of the system. The complete equations of motion can be 
written as a variation about the equilibrium point (x 0 ,0,/o) 


u = x - Xq 

v = jc (11) 

W = I -Iq 


where 


d_ 

di 


< 


u 

v 

w 


) 



A = 


0 1 

T ff 0 

0 2 

|_ a 




( 12 ) 


(13) 


where 



Cl c ’]0 


a = /j.qN 2 A 


(14) 

(15) 


Note, the equilibrium values ( xq,Iq,Eq ) are related and consistent values 
a fixed point. 



must be used so that the position is actually 


(16) 
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(17) 


/ - Eo 

Io ~T 

The various terms in Eq (13) can be identified as representing the magnetic force, the voltage drop due to back 
emf and the voltage drop due to inductance. 

The controllability of this system can be determined from the rank of the controllability matrix Q. 


Q — [B AB A 2 B ] 

Substituting from Eq. (13) results in a matrix that has full rank. The total dynamic matrix turns out to be 


(18) 
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(19) 


The system is stable if 


< R 


> 2R 

A£g ( 

a 

< 2 R 

| 

a 


The first constraint can be interpreted as requiring the coil to behave with positive resistance. The minimum value of 
Ci is that required to overcome the effective negative stiffness of the magnetic force. The third condition enforces an 
upper bound on Ci which is dependent on Cz. This implies that the system cannot be stable without some velocity 
feedback. Apparently, the damping induced by the back emf is cancelled by the inductive lag of the coil, leaving an 
unstable system unless velocity feedback is used to add dissipation. 

Given values of Ci,C2,C3, the system may be linearized about the other equilibrium point ( the one at lower 
gap). This point is found still to be unstable. Fig. 3 shows a sketch of the response in 3-D state space. The 
stabilized equilibrium point is a focus in two dimensions and has a stable subspace in the other dimension. The 
unstable equilibrium point is a saddle point in two dimensions and has a stable subspace in the third dimension. The 
presence of the unstable equilibrium point close to the stable one significantly affects the domain of convergence. 

Rather than developing this problem ftirther ( including for example, the effect of eddy currents) discussion will 
be shifted to a two magnet configuration. 

TWO MAGNET VERTICAL CASE 

Consider a bearing consisting of only the top and bottom magnets in Fig. 1, with x constrained to be zero. 
This problem will study how the eddy currents affect the system, and so Appendix B will be used to represent the 
magnetic forces. The rotor is suspended by the difference in magnetic forces due to the top and bottom magnets. 
The stationary steady state voltages are such that the gaps are equal top and bottom, and the magnetic flux density 
is well below the saturation value of 1.5 - 2.0 webers/m 2 . 

Since the system is constrained in 
The governing differential equations 


The 8 state vector reduces to (y,y,/i,/ 3 ). ZF mag becomes F magl -F mag3 . 


x, the drag force due to eddy currents and the resulting coupling is neglected 
are 
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E 03 4- <5^3 —^3^3 
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2y/i cos 0i 


( 21 ) 


2(A — y cos0i) +P\ 
2 y/ 3 cos 03 
2(A — y cos 03 ) +/J3 
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The voltages £oi,Eo 2 have been adjusted so that the equilibrium point is y = 0 with no spin speed. It then 
follows that 

, _ £ oi 

701 -7FT 


^03 ~ 


Ep 3 
R 3 


( 22 ) 


£ R m agl Rm ag3 

A state variable feedback controller is designed for this case using the linearized set of equations at zero speed 
(£2 = 0.0). Because eddy currents are not a factor at zero speed , the formulae in Appendix A can be used ( giving 
a significant simplification over Appendix B ) . The linearized system can be written as 

tmAZ+, (to,) <23) 

A is the Jacobian, B is the control matrix containing the inverse of the inductances of the two circuits at equilibrium 
and ( SEi SE i) T is the control vector for the perturbation voltges to the magnets. 

This is a two input - single output system and unlike the previous single magnet case, there are many gain 
matrices K yielding the same eigenvalues. The solution to this control problem is discussed in [18]. However, the result 
is a gain matrix K such that the poles of the closed loop system ( at zero speed ) are placed at -1000, -1000, -100± 
j 100. 



The complete closed loop system can now be considered with speed as a parameter. The equilibrium points 
are found by the numerical solution of the right hand side of Eq. (21) with ( y = 0). It turns out that there are 3 
equilibrium points. This can be understood from a graph of the net magnetic fource versus y as in Fig. 4. The point 

at y = 0 is a stable equilibrium point. However, points at y ~ 4.4 10“ 4 and y 8.4 10 -4 are also equilibrium 

points. 

Furthermore, as speed increases, the force curve changes. The central equilibrium point drifts downwards, and 
the slope there decreases. Eventually, at sufficiently high speed, the central and lower equilibrium points coalesc and 
disappear. This is shown in Fig. 5. 

The catastrophe occurs at £2 = 47000rpm for the values of parameters in this paper. 

However, before that occurs, a more interesting dynamic phenomenon is observed. The dynamic equations 
can be linearized at each (shifting ) equilibrium point. Since the equilibrium points must be found numerically, the 
linearization was also performed numerically using a central difference scheme in each state variable. The eigenvalues 
of the resulting dynamic matrix are shown in Fig. 6 as a function of the shaft speed £2. 

At zero speed, the eigenvalues are very near the design points -1000, -1000, -100 ± j 100. The difference is 
because magnetic circuit theory ( Appendix A) is used in the zero speed controller design, but the full eddy current 
field theory solution ( Appendix B ) is used when linearizing and calculating eigenvalues. The system is stable ( 
in a small neighborhood of the central equilibrium point ) up to a value of £2 = 4G406rpm, at which time a Hopf 
bifurcation (flutter instability) occurs. 

The system loses stability through a single pair of eigenvalues crossing the imaginary asis at £2 = 40406rpm. In 
such a case, a limit cycle must grow from the bifurcation point. This can be either a stable limit cycle which increases 
in amplitude for £2 > 40406rpm or an unstable limit cycle which closes down around the focus as £2 -» 40406rpm. 
A standard bifurcation analysis was performed using the program BIFOR2 [19]. A stable supercritical limit cycle is 
predicted. The amplitude is shown in Fig. 7 along with the limit cycle actually found by simulation. It should be 
noted that the actual phase space is four dimensional and only a 2-D subspace is shown. 

The reader is cautioned that the study was not extended to determine what ( if anything) happens to the limit 
cycle when the two equilibrium points coalesc and disappear at 47000 rpm. 


FOUR MAGNET CASE 

Each coil is excited by initial constant voltages denoted by E oi , i = 1,2, 3, 4, and the steady state magnetic forces 
(at £2 = 0) are such that 

y^Fmag = -mghy. (25) 
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The equilibrium points for this system at any spinning speed can be determined by equating the right hand side of 
Eqs. (5) to zero. That is, 

0 = /(X,„Q) (26) 


give the equilibrium points of the system. Eq. (26) represents a set of coupled nonlinear algebraic equations which 
has to be solved by a numerical technique. But, at zero speed, the equations decouple in x and y giving two separate 
problems, one of which has been discussed briefly. Without feedback control, this system has a single equilibrium 
point. It is possible to track the movement of this equilibrium point as speed increases [18]. 

To find the nature of the equilibrium point at any speed Q, the differential equations are linearized and studied 
for their behavior in the neighborhood of the equilibrium point. As expected, the eigenvalues of the linearized system 
show that the four magnet bearing without active controls is inherently unstable. 

A state variable feedback controller is added to the system to stabilize the bearing and provide it with suitable 
damping and dynamic characteristics. The effect of a gain matrix which is speed dependent is beyond the scope 
of this paper, and the control system is designed using the system characteristics at Q = 0.0. At zero speed, the 
equations of motion are uncoupled, and two 4th order control problems can be solved rather than one 8th order 
problem. State variable feedback controllers were designed for each subsystem to place the closed loop poles at -1000, 
-1000, «100±jl00. This placement provides dominantly second order response over a frequency range of 0 to 5000 
rpm. 


Note that because the system decouples, all states are not fed back to the control inputs for each magnet. For 
example, the variations of the gap length, velocity, and currents of the horizontal magnets are not being fed back to 
the vertical magnets and vice versa. Another advantage is that designing a state variable feedback controller for a 
four-input two-output system is more involved in terms of computation that two versions of a two-input single-output 
system. The feedback gain matrix K for this four magnet bearing is given in Table 3. 
where 

r se ,i 


se 2 

5E 3 

. se 4 . 


= [K-JX 


(27) 


It is not possible to derive analytical expressions for the equilibrium points of the system. However, the pattern 
of 3 equilibrium points for a pair of magnets in one dimension is extrapolated to an expectation of 3 2 equilibrium 
points in 2 dimensions. This has been confirmed by numerically solving the equations. As in the simpler cases, in 
designing the control system, the choice of pole assignments for the central equilibrium point (and resulting feedback 
gains and stiffness) affects the location of the surrounding equilibrium points. A faster, stiffer system (poles further 
into the left half of the s-plane) pulls the outlying equilibrium points closer together, as shown in Fig. (8). 

As the spin speed increases, the equilibrium points shift due to changes in the magnetic forces, and in this case 
including the coupling due to eddy current drag forces. Furthermore, because the dynamics are affected by the gap 
lengths, the dynamic coefficients have changed, and the stability of the system is affected. Fig. 9 shows how all 9 
points move within the clearance space. At low speeds not much happens but above 10,000 rpm the points begin to 
shift. At about 15,000 rpm points 3 and 9 coalesce and disappear. 

At any running speed, the system can be linearized about any of the equilibrium points. Since the equations 
are rather unmanagable, the linearization was performed using a central difference method to find the derivatives of 
each function with respect to each state variable. The result is 


Z = AZ. 


(28) 


where the Z column vector contains perturbations about the speed dependent steady-state vector X. Eigenvalues of 
A give stability information about each point. Points 2-9 ( the outlying points ) were always unstable. 

The eigenvalues corresponding to the central equilibrium point ( 1) as a function of the spinning speed are plotted 
in Fig. 10. The system is stable in a small neighborhood around the equilibrium point up to about 9306 rpm. As 
the spinning speed is further increased, a pair of complex conjugate eigenvalues crosses the imaginary axis with all 
the other eigenvalues still in the left half of the complex eigenvalue plane. This means that the system undergoes a 
classical Hopf bifurcation to flutter instability at the critical speed of 9306 rpm. 

Because of the existance of the outlying unstable equilibrium points, a supercritical limit cycle was expected. 
However, a conventional Hopf bifurcation analysis on the full set of nonlinear equations, using BIFOR2, indicated 
that an unstable limit cycle will enclose the central equilibirum point for speeds less that Qc r . Fig. 11 shows three 
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transient responses at G = 2500 rpm, calculated by numerical integration of the complete set of nonlinear 
equations. Two converge to the equilibrium point and one diverges. The domain of convergence was approximated 
by slowly incrementing .x(O) and y(0) and noting convergence or divergence. The domain of convergence is marked 
in Figs. 11,12. This is an overly simple approach because changes in the initial conditions of the remaining 6 state 
variables are not explored. However, it does lend credence to the prediction of an enclosing unstable limit cycle. 

Forced response of this system is discussed in detail in [18]. The global resopnse is quite different from what 
might be expected from a linear analysis. The linear eigenvalues at ft = 2500 rpm indicates a damping ratio of 

0.7. However, the orbit resulting from a periodic excitation of 25N only converges quite slowly (damping ratio < 
0.05). Further discussion is beyond the scope of this paper. The topic is opened here only to point out the pitfalls 
in extrapolating a linearized analysis. 

CONCLUSION 

Without considering eddy current effect, there is no speed dependent term in the constitutive model for a 
magnetic bearing. Eddy current cause a loss of effective lift which, viewed as as external load, causes a classical sag 
for a proportional controller. More importantly, eddy current drag causes coupling between (x,y). Any change in 
the x gap affects the drag force in the y direction. Straightforward proportional control produces extra equilibrium 
points as the applied voltage is driven negative but the force remains attractive. These points are unstable but affect 
the global response of the system. In fact, designing a stiff system draws these points quite close to the equilibrium 
point. The equilibrium point is surrounded by an unstable limit cycle. Furthermore, as speed increases, the unstable 
limit cycle shrinks until the central point loses stability in a classical Hopf bifurcation. Above critical velocity, there 
is no equilibrium response. 

The loss of lift is about 30% at 12500 rpm, which agrees approximately with experimental results of Yamamura 
[20]. The expression for the drag force under a single sinusoidal field density wave is the same as that obtained by 
Meisenholder [1]. 

The system may also lose stability by a simple catastrophe, the coalescing and joint annihilation of two equi- 
librium points. For the parameters used in this paper, this happens after the Hopf bifurcation. However, further 
parameter studies are needed to determine if this might occure before the Hopf bifurcation. 

Further study is needed in designing state variable controllers. The zero speed design point is useful because it 
decouples the system, reducing the order. However, many other types of control can be envisioned. 

The size of the enclosing unstable orbit is quite small, as shown in Fig. 12. A transient excitation could rather 
easily throw the shaft outside this into an unstable response. Also, a periodic excitation such as rotating unbalance, 
produces an orbit about the central equilibrium point. Simulation has indicated that if this orbit is large enough 
to touch or exceed the enclosing unstable limit cycle, the system becomes unstable. The loss of stab iliity is quite 
complicated, with bifurcations of the periodic orbit and possibly bifurcation to chaos. Discussion of this is beyond 
the scope of this paper. 

Further problems which should be investigated are the effects of flux saturation in the magnetic material. Also, 
the model developed in Appendix B is applicable to repulsion type electromagnets, which would run at high speed 
and likely require superconductivity to implement. Stability of these bearings has not yet been addressed. 

A final comment is the topic of bearing coefficients ( equivalent stiffness and damping matrices ). The situation is 
different from that of journal bearings, it is not possible to determine the 8 classical coefficients just by differentiating 
the force expressions with respect to x,y,x,y. The dynamics of the electrical components must also be incorporated. 
However, in all cases, the response in one ’mode' was very fast( Fig. 3) indicating that some subspace reduction 
might be possible. 


APPENDIX A 

ESTIMATION OF MAGNETIC LIFT FORCE 

Magnetic circuit theory can be used to approximate the magnetic lift force of a single magnet, but 
eddy currents must be neglected. The following assumptions are made in deriving the expressions for the 
magnetic lift force: 

1. Field fringing is neglected. 

2. Magnetization curve is linear (B = /i//). 

3. Magnetic flux density B and field intensity H are uniform over cross-sections of the core, gap, or mass. 
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An electromagnetic circuit is considered whose elements are gap, core, and suspended mass, as shown in Fig. 
A.l. Each element has constant crossection A,- and length L*. The magnetic flux (p is assumed constant throughout 
the circuit, and T is the total magnetomotive force within the circuit elements. 

The density relationship <p = BA and the constitutive law B = /iH can be used to express field intensity in the 
ferromagnetic material in terms of the field intensity within the air gap. 

By definition r= / Hdl Integrating around the circuit and equating T to the current linked (NI) results in an 
equation for field intensity within the air gap. The field energy is determined within each element w = \jxH and the 
total field energy is obtained by a volumn integral over all the elements. 

By definition, the force is the rate of change of stored field energy with respect to the mechanical displacement. 
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In addition, the total flux <p, the magnetic field density B, the magnetic flux linkage X, and the inductance L can be 
expressed in terms of the gap length x and current I as: 
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APPENDIX B 

ESTIMATION OF MAGNETIC FORCES 
INCLUDING EDDY CURRENT EFFECTS 

Motion of a conducting material through a magnetic field will cause eddy currents to be generated within the 
material. These eddy currents will produce an additional magnetic field and change the net force acting across the 
air gap. The complete eddy current analysis for 4 magnets as shown in Fig. 1 is analytically intractable. (It is 
probably amenable to finite element techniques.) An approximate solution is developed in Ref [18], and only the 
assumptions, general method, and results will be presented here. 

As the shaft starts spinning, the eddy currents tend to repel the applied magnetic field and the skin depth of 
penetration becomes very small. This motivates a semi infinite assumption in the radial direction. 

First, the problem will be unwrapped and considered as periodic on a half-space. However, before net forces are 
calculated, the surface tractions predicted by Maxwell’s stress tensor will be wrapped around a circular shaft. 

Second, each magnet will be considered separately, and the magnetic field for each magnet can be determined 
individually. The net force of each magnet is then determined, leading to 4 vector forces which are then summed 
vectorially. An alternative (more complex) solution is developed in Ref. [18] to find net magnetic field for all four 
magnets as a single system (simultaneously). The net magnetic field of all 4 magnets may be determined, and a 
single force predicted. However, this approach requires the assumption that all the gaps are equal. At low spinning 
speeds there is no difference between the two methods. Only at very high speeds do the two differ. (For parameters 
in this paper, 4% at 10 5 rpm). The simpler technique has been used here. 

The square wave applied flux density is expanded as a Fourier series 



or 

B(x) - ^BiSin(kix) (B.2) 

i 

The field density and current density distributions within the moving material solve a linear problem, and hence 
the principle of superposition can be invoked and each harmonic handled separately. After finding the field density 
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inside the whole region, the forces acting on the rotor can be determined by calculating the Maxwell’s stress tensor 
and integrating it over the surface area. 

The following assumptions are made to simplify the analysis: 

1. The conductive plate is infinite in x,z and positive y directions. 

2. The conductivity o and permeability fi are constants. 

3. The field problem is two dimensional. 


The equation that describes the distribution of the magnetic field in the conducting medium is derived from 
Maxwell’s field equations [7]. 

1 ■> dB -» -> 

V 2 B+~. r- = Vx(Vx£) (B.3) 

jucr dt 

where B = B x e i + B y e 2 (ei and e 2 being unit vectors in x and y directions respectively). The y component B y 
is determined from the y component of (B.3) The remaining component B x can be determined from the relation 
V B =0. The magnetic field is driven by the applied magnetic field density, and so solutions with the same traveling 
wave dependence on (x,t) are assumed. That is, it is assumed that the flux density takes the form 

B = [B x {y)e, + B y (y)e 2 ¥ <**""'> (B.4) 

The solution form in the y-direction is where 
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Eqn (B.5) can be used to predict the skin depth. 

The solution domain is divided into two regions, denoted by subscripts 1 and 2 respectively: 

Region (1), the air gap where a = 0(0 < y < A ); and 
Region (2), the moving conducting medium (A < y < oo). 

The solutions within each region have two constants of integration, determined from the following 4 boundary 
conditions: 


(1) one is set to zero because the solution cannot grow as y -^oo. 

(2) the applied flux density at y = 0 

(3) At the interface y = A , the condition of conservation of flux V B = 0 is invoked. Using the divergence 
theorem, this leads toB y i = B y2 . 

(4) at the interface, V xH = 0. Using Stokes theorem, this leads to H x \ - H x2 . 

Hence, the flux density distribution throughout both regions can be determined. 

The forces acting on the conducting medium are calculated by Maxwell’s stress tensor [8]. For magnetic problems 
with currents and no charges, the forces acting on a body are given by 


F — f -[BB h — }-B 2 h]dA (B.6) 

J E ll 2 

where Z is any closed surface surrounding the body and not containing any other body and B is the value of the 
field on the closed surface., Choosing a closed surface Z such that it extends from -oo tooo and includes only the 
conducting medium, the integration is carried out with h = -e 2 . 

The complete flux density distribution in the whole region of Fig. B.l due to all the applied sinusoidal waves is 
determined by superposing the individual fields. Each component of the field (B x ,B y ) is an infinite series in sine or 
cosine terms. The value of the B field at the interface, which is required for calculating the forces, is calculated by 
substituting y = A and letting B x = B x \ and B y = B y \ . Substituting for B and evaluating the integrals, the forces 
per unit area acting on the material turn out to be 
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Eqs. (B.7) and (B.8) are the expressions for the drag and lift force per unit area acting on the moving conducting 
medium. As expected, when the currents are not induced in the slab (V=0) (S— 0), there is no drag force (F x ) and 
Fy is the magnetic attractive force. There is an optimum value of V at which the maximum force per unit area F x 
is produced. The lift force decreases as V is increased, and at some value of V the force becomes zero, and at high 
values acts in the opposite direction (repulsion). At very high values of V, there is no drag because all the flux is 
excluded from the material and the repulsion force reaches an asymptotic value irrespective of the permeability of 
the material. 

It is interesting to note that the integral and the summation are interchangable in order. That is, the force for 
one component of the fourier series can be determined and then summed or, as in the preceding, the fields summed 
and the force determined. This is perhaps surprising because the problem is nonlinear, but the infinite series for B x 
and By are made up of sine and cosine terms which are orthogonal to one another and all the cross terms drop out 
during integration over one period. 

The total flux can also be calculated and compared with that predicted by magnetic circuit theory (Appendix 
A). The more detailed solution is about 8% lower, showing the effects of magnetic circuit assumptions (uniform field 
density and no leakage in air gap). 

For the magnetic bearing, the forces acting on a rotating shaft are calculated by wrapping one period of the B 
field distribution back onto the circular shaft. Choosing a closed surface X on the circumference of the rotor and 
simplifying the integral in Eq. (B.6) give the forces acting on the rotor per unit width as 
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By = By\(y ) = i€ ^ k X 
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(B.14) 


The effect on the values of the forces of the number of terms used in the fourier series was investigated. The 
magnetic forces were calculated at different spinning speeds using 10, 25, 50, 100, and 500 terms in the series. There 
is very little change in the results when the number of terms used in the series is 50 or more. There is approximately 
3% change in the lift force when the number of terms is increased from 10 to 25 or from 25 to 50. The change is 
much less in the case of drag force calculations. Striking a compromise between these two, and to save computer 
time, 25 terms were used in all the calculations in this paper. 
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TABLE 1 

Single Magnet Suspension Parameters 


Mass M 

Acceleration due to gravity g 
Area of cross section A 
Initial voltage E 0 
Resistance of coil R 
Path length in coreL] 

Path length in Mass L 2 
Permeability of free space /r 0 
Relative permeability of core material p r \ 
Relative permeability of mass material p r z 
Number of turns N 


= 8.91 kg 
= 9.80665 m/sec 2 
= 0.00025 m 2 
= 5.0 volts 
= 3.0 ohms 
= 0.28 m 
= 0.05 m 
= 4 kE — 07 H /m 
= 10000.0 
= 100.0 
= 800 


FEEDBACK GAIN MATRIX: Poles at -1000 , -100 ± jlOO 
K = [261013.13 1486.86 -92.44] 


409 



TABLE 2 

Two Magnet Bearing Parameters 


ROTOR 

Diameter of rotor D 

Mass of rotor M 

Path length in rotor L 2 

Relative permeability of mass material 

Electrical conductivity of mass material a 


MAGNET 

Pole angle 

Path length in coreLi 

Relative permeability of core material fi r i 

Area of pole face A 

Width of pole face W 

Number of turns N 

Distance of pole corner a 

Distance of pole corner b 

Initial voltage (top) £ 0 i 

Initial voltage (right) E o 2 

Initial voltage (bottom) E 03 

Initial voltage (left) £04 

Resistance of coil R 


= 0.15 m 
= 15.6 kg 
= 0.05 m 
= 100.0 

= 1.0 E+07 mho/m 


= 15 degrees 
= 0.28 m 
= 10000.0 
= 0.00025 m 2 
= 0.008333 m 
= 800.0 
= 0.0075 m 
= 0.0375 m 
= 7.5 volts 
= 2.5 volts 
= 2.5 volts 
= 2.5 volts 
= 3.0 ohms 


TWO MAGNET FEEDBACK GAIN MATRIX 
Poles at -1000, -1000, -100+ jlOO 


' 4570751.1 11052.9 -565.5 -2200.5' 

1523583.7 3684.3 -188.5 -733.5 . 


TABLE 3 

Feedback Gain Matrix for Four Magnet Bearing 


4570751.1 

0.0 

1523583.7 

0.0 

11052.9 

0.0 

3684.3 

0.0 

0.0 

7185390.6 

0.0 

7185390.6 

0.0 

18021.2 

0.0 

18021.2 

-5675.5 

0.0 

-188.5 

0.0 

0.0 

-3500.0 

0.0 

-3500.0 

-2200.5 

0.0 

-733.5 

0.0 

0.0 

-3677.7 

0.0 

-3677.7 
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Figure 5: Locus of Equilibrium Points, with Feedback, as Speed Varies 



Figure 6: Root Locus for Two Magnet Bearing, as Speed Varies 
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Figure 7: Limit Cycle Comparison for Two Magnet Bearing 
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Figure 8: Equilibrium Points at £2 = 0 for Three Different Pole Assignments 
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Figure 9: Locus of Equilibrium Points as a function of Speed 
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Figure 10: Root Locus as Speed Varies 
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Figure A.l: Ferromagnetic Rotor under a Single Electromagnet 



Figure B.l: Semi-Infinite Plane under a Series of Electromagnets 
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DRAG FORCE ( N ) 


EOOY CURRENT EFFECTS 



Figure B.2: Lift Force as a function of Spin Speed 



Figure B.3: Drag Force as a function of Spin Speed 
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